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Purpose: We analyze a recently proposed polyenergetic version of the simultaneous algebraic 
reconstruction technique (SART). This algorithm, denoted pSART, replaces the monoenergetic 
forward projection operation used by SART with a post-log, polyenergetic forward projection, while 
leaving the rest of the algorithm unchanged. While the proposed algorithm provides good results 
empirically, convergence of the algorithm was not established mathematically in the original paper. 

Methods: We analyze pSART as a nonlinear fixed point iteration by explicitly computing the 
Jacobian of the iteration. A necessary condition for convergence is that the spectral radius of the 
Jacobian, evaluated at the fixed point, is less than one. A short proof of convergence for SART is 
also provided as a basis for comparison. 

Results: We show that the pSART algorithm is not guaranteed to converge, in general. The 
Jacobian of the iteration depends on several factors, including the system matrix and how one 
models the energy dependence of the linear attenuation coefficient. We provide a simple numerical 
example that shows that the spectral radius of the Jacobian matrix is not guaranteed to be less than 
one. A second set of numerical experiments using realistic CT system matrices, however, indicates 
that conditions for convergence are likely to be satisfied in practice. 

Conclusion: Although pSART is not mathematically guaranteed to converge, our numerical 
experiments indicate that it will tend to converge at roughly the same rate as SART for system 
matrices of the type encountered in CT imaging. Thus we conclude that the algorithm is still a 
useful method for reconstruction of polyenergetic CT data. 


I. INTRODUCTION 


The algebraic reconstruction technique (ART) or 
Kaczmarz’ method [1, 2] is a well-known method for ap¬ 
proximately solving systems of linear equations, 

Ax = b, (1) 


where x is a column vector of size n, b is a column vector 
of size to, and A is to x n. The method has a long his¬ 
tory in CT image reconstruction, in which b represents 
the post-log measured projection data, x represents the 
image to be reconstructed, and the (i,j )th entry of A 
represents the length or area of intersection of the ith 
ray with the jth pixel of the image. The values in x rep¬ 
resent the averaged linear attenuation coefficient (LAC) 
within each pixel, and are usually constrained to be pos¬ 
itive. Since the model is linear, it is implicitly assumed 
that the X-ray beam used to generate the data is mo¬ 
noenergetic, and the reconstructed values correspond to 
the LAC of the tissue at that energy. 

Beginning from an initial estimate x^°\ ART itera¬ 
tively projects the current iterate, onto the hyper¬ 
plane defined by one of the m equations defined by (1). 
The ART iteration is given by 


x (fc+i) = x (fc) _ 


c (fc) ) - bi 


(a,, a*) 


( 2 ) 


where a* is the ith row of A, (•, •) denotes the dot prod¬ 
uct, and i is chosen to be (k mod to) + 1 in the classical 
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version of the algorithm. Some references consider this 
step to be a sub-iteration, and a full iteration of ART 
to consist of successively applying this iteration for all 
m equations. If the system is consistent, the iteration 
is guaranteed to converge to a solution, while in the in¬ 
consistent case the sub-iterations converge to a limit cy¬ 
cle [3], 

The simultaneous ART (SART) method [4, 5] is a vari¬ 
ant of ART in which the corrections generated by the 
ART sub-iterations (2) are combined and applied simul¬ 
taneously. The iteration can be expressed concisely in 
terms of matrix operations as 

x (fe+1 ) = x (fc) - DA t M (Ax (fe) - b) , (3) 

where D and M are diagonal matrices, with 


Djj „ , f3j — |Ofcj |,j — 1 ... i 


M u = 


ft 

1 

li' 


k =l 


li = y) \a,ik\,i = 1.. .to. 


(4) 


fc=i 


In other words, /3j is the 1-norm of the jth column of A, 
and 7 j is the 1-norm of the ith row. 

The primary advantage of SART over ART is that 
SART is less sensitive to noisy data [4]. It should be 
noted that the iteration (3) computes an update using 
projection data corresponding to all views simultane¬ 
ously, while the original SART algorithm [4] computes a 
sequence of updates using projection data corresponding 
to only one view of the object at one time, to accelerate 
convergence. For simplicity of analysis we will only con¬ 
sider equation (3), as this was the used as the basis for 
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the polyenergetic algorithm that we will analyze. Like 
ART, SART has been proven to converge [6, 7]; see also 
the Appendix to this paper. 

As mentioned, both ART and SART implicitly assume 
that the data are generated from a monoenergetic X-ray 
beam, i.e. that 


Pi = /exp(-(a l ,x)), (5) 

where Pi is the measured intensity of the ith beam, and 
/ is the blank scan intensity (assumed to be independent 
of i). Taking the log of the data and rearranging terms 
then gives a linear system equivalent to (1). 

-hiy = (0i,x). (6) 

In practice, however, the X-rays generated by clinical CT 
hardware are usually polyenergetic. A typical model for 
polyenergetic X-ray measurements is 

Pi = /"/(e) exp (-(a*, n(e))) ds, (7) 

where e refers to the energy of an incident x-ray, /(e) is 
the initial intensity of the beam corresponding to that 
energy (i.e. the beam’s spectrum) and fi(e) is a vec¬ 
tor of attenuation coefficients, whose values depend on 
X-ray energy. This system of equations can no longer 
be linearized, and it is well-known that reconstructing 
an image from this data using conventional means (such 
as ART, SART, or filtered back projection) produces an 
image containing beam hardening artifacts [8]. This mo¬ 
tivates the need for polyenergetic iterative reconstruction 
algorithms, e.g. [9-18]. 

The polyenergetic SART (pSART) algorithm [18] was 
recently proposed as one such reconstruction technique. 
In this approach, the beam spectrum is discretized into 
h energy levels Eh, with weighting terms Ih computed 
to approximate the continuous spectrum. The vector¬ 
valued function /x(e) is then modeled as a function of a 
vector t, representing the attenuation map of the object 
at a reference energy, Eo, which was chosen to be 70 keV. 
This function makes use of tabulated energy-dependent 
LAC values for some suitable reference materials, such as 
air, fat, breast, soft tissue and bone. Letting tj denote 
the LAC value for pixel j at the reference energy, the 
LAC of that pixel for all other energies e is then given 
by 


soft tissue, with the weighting determined by the values 
at the reference energy. 

One can then define a polyenergetic forward projection 
operator, V : R n —¥ R m , which acts on t: 


[P(t)]i = 4 exp (-(a,, p(t, E h ))) . (9) 

h 


The pSART iteration is then defined as 


t ( fc +i) = t ( fc )-DA T M(- ln -p(V fe) ) +ln(p)), (10) 


with D and M defined as in (4). Equations (9) and (10) 
are equivalent to equations (14) and (15) in Ref. 18, al¬ 
though our notation is somewhat different. The only 
difference from the SART iteration (3) is that the log of 
the monoenergetic forward projection has been replaced 
by the log of the polyenergetic forward projection. The 
algorithm produces a single attenuation map of the ob¬ 
ject, t, with LAC values corresponding to the reference 
energy. 

In Ref. 18 the authors state that this modification 
solves the problem of inconsistency between the polyen¬ 
ergetic data and the monoenergetic model implicitly as¬ 
sumed by the conventional SART approach. While it is 
clear that a vector t satisfying p = V(t) is a fixed point 
of this iteration, this does not guarantee convergence of 
the algorithm. Experimental results indicated that the 
method is effective, however. In the next section we ana¬ 
lyze the convergence of pSART as a fixed point iteration. 


II. CONVERGENCE OF PSART 

We first consider the SART iteration (3). The iteration 
can be written in the form 

x (fc+ i) = T yS k) + C , (11) 

where T = I — DA T MA and c = DA T Mb. One can 
show that the spectral radius of T (the magnitude of 
its largest eigenvalue), denoted p(T), is strictly less than 
1 (see Appendix). This guarantees convergence of the 
algorithm to a solution, if the system is consistent. It 
has been proven that if no exact solution exists, SART 
converges to a weighted least-squares solution [6, 7]. 

The pSART algorithm is a nonlinear fixed point itera¬ 
tion. We write the iteration as 


( , x _ K--t-i(go) ~ tj]p, k (e) + [tj - p k (e o)W+i(e) 
m ' } p k+1 (E 0 )-Me o) 

( 8 ) 

where /ifc(e) and pk+ii^) are the tabulated, energy- 
dependent LAC functions for the two base materials with 
LAC values adjacent to tj at the reference energy. So for 
instance, if the value of tj is between the LAC for soft 
tissue and the LAC for bone at the reference energy, then 
its LAC at all other energies is obtained by linear inter¬ 
polation between the corresponding values for bone and 


t ( fe ) _ DA T Mf (t (fe) ) , 

(12) 

F (t (fc) ) 

(13) 


where /(t) = — ln['P(t)] + ln(p). Assuming that a solu¬ 
tion to the system of nonlinear equations exists (which 
corresponds to a fixed point of the iteration), a necessary 
condition for convergence is that the spectral radius of 
the Jacobian matrix of F, p{Jp), must be less than one 
when evaluated at the solution. Note that this does not 
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guarantee that the algorithm converges to the solution 
from any starting point; simply that it cannot converge 
if the condition does not hold. Let t* denote a solution 
to the nonlinear system of equations p, = (P(t)]j. A 
straightforward calculation gives 

J F (t*) = I - DA T MJ f (t*), (14) 

where the (i,j)th element of the m x n Jacobian matrix 
J/{ t) is given by 


dfj 

dtj 


(t) 


C^ij 

WWi 


I h exp (—(oj, fi(t, e h ))) 

. h 


It follows from (8) that 


dp 

~dt 


{tj , Sh ) 


(15) 


Mk+i(<0 -Mfc(e) 

o, f \ i \ ? (lu) 

dt £ifc+i(eo) — Mfe( £ o) 

where pu+i and pk again refer to the base material LACs 
that are adjacent to t at the reference energy e 0 . This 
implies that for a fixed value of e, this partial derivative 
is piecewise constant, with discontinuities when t is equal 
to one of the reference LACs at £q. 

A general analysis of the spectral radius of Jp{ t*) is 
difficult as it has a complicated dependence on the system 
matrix A, the spectrum Ih, and the choice of base ma¬ 
terials. We now show with a numerical experiment that 
one cannot guarantee that p(Jf{ t*)) < 1, in general. 

We consider a simple example where m = n = 2, illus¬ 
trated in Fig. 1. The object consists of two pixels of size 
lxl cm, with LACs of ti = 0.1 cm -1 and t 2 = 0.16 cm -1 
at a reference energy of 70 keV. We first consider the case 
of a monoenergetic 70 keV beam with intensity 7=1. 
The first beam travels through both pixels horizontally, 
while the second beam has a length of intersection of 
roughly 0.28 cm with the first pixel and 1.13 cm with 
the second pixel. After taking logarithms, we obtain the 
2x2 linear system of equations 


1 

1 


tl 


'0.260' 

0.28 

1.13 


^2 


0.209 


which can be solved using ART or SART. Fig. 2 shows 
the progression of both ART and SART for this system, 
starting from an initial guess of zero. One can see that 
ART sequentially projects onto the two lines defining the 
equations, while the iterates generated by SART follow 
a path between the two equations. 

We now consider the case of a polyenergetic x-ray spec¬ 
trum. The LACs of the two pixels at the reference en¬ 
ergy determine the LACs at all other energies e according 
to (8). The discrete spectrum used in this experiment 
consists of eleven energy bins, obtained from a continu¬ 
ous 130 kVp spectrum generated using the Siemens Spek- 
trum online tool [19, 20]. This spectrum and the atten¬ 
uation curves for the two materials are shown in Fig. 3. 
The attenuation values for the reference materials were 


P 2 = exp(—0.209) 

71 



/ 

/ 

/ 

/ 

/ 

/ 

h = 0.1 / 
/ 

/ 

/ 

t 2 = 0.16 


pi = exp(—0.260) 

> 


FIG. 1. Example problem consisting of two pixels. 


Monoenergetic case 





FIG. 2. Convergence of ART & SART for the monoener¬ 
getic experiment. The black lines correspond to the system 
of equations (17), red line correspond to the sequence of iter¬ 
ates generated by ART, blue line to the sequence of iterates 
generated by SART. 


obtained from Ref. 21. Under the polyenergetic model, 
we obtain a nonlinear system of equations: 

Pi = ^2lhexp(-p(h,£ h ) - p{t 2 ,e h )) 
h 

ss exp(—0.314) 

p 2 = ^I h exp (—0.28/i(ti, e h ) - 1.13 p(t 2 , £ h )) 

h 

« exp(—0.253) 

The two beams undergo more attenuation than in the 
monoenergetic experiment (17) because the spectrum 
contains a higher proportion of X -rays with energies less 
than 70 keV. In Fig. 4, this nonlinear system of equa¬ 
tions is illustrated with black curves. As either t\ or 
t 2 increases, the LAC value of the pixel at lower ener¬ 
gies increases rapidly, meaning that the coefficient in the 
other pixel must decrease rapidly to compensate. Thus 
the curves are bent. For the sake of illustration, we have 
implemented a polyenergetic version of ART (denoted 
pART) that is analogous to pSART. The red and blue 
lines show the progression of the pART and pSART it¬ 
erations, respectively. For this experiment (left figure), 
one can see that both iterations converge to the solution 
of the nonlinear system of equations. 

We now give a case where the iteration fails to con¬ 
verge. In the right figure of Fig. 4, the progress of the 
two iterations is shown for a slightly modified experiment, 
where t 2 was changed from 0.16 cm -1 to 0.24 cm -1 . All 



















4 



Spectrum (normalized) 



FIG. 3. X-ray spectrum and attenuation curves used for the 
polyenergetic experiment. Left: Continuous spectrum (blue 
line) and discrete energies (red crosses) used for the summa¬ 
tion. The spectrum has been normalized to have an integral 
of 1. Right: Attenuation curves for the base materials as well 
as the interpolated curves for values of t = 0.1, 0.16 and 0.24 
(red dashed lines). The reference energy of 70 keV is indicated 
by the dashed black line. 


Polyenergetic case 


Polyenergetic case 




FIG. 4. Convergence of pART and pSART for two polyener¬ 
getic experiments. The black curves correspond to the nonlin¬ 
ear system of equations representing the polyenergetic mode, 
while the red and blue lines correspond to the sequence of 
iterates generated by the pART and pSART algorithms, re¬ 
spectively. In the first experiment (left plot) the LAC values 
at the reference energy are ti = 0.1 and t 2 = 0.16, while in 
the left plot, t 2 = 0.24. 


other parameters of the experiment were the same as 
before. It is apparent that both iterations fail to con¬ 
verge to the solution; pART appears to exhibit chaotic 
behaviour about the solution, while pSART converges to 
a two-cycle. Even if the iteration is started very close to 
the solution t* = (0.1,0.24) T , both pART and pSART 
diverge. 

A direct computation of the 2x2 Jacobian matrix (14) 
for these two experiments reveals that the spectral radius 
of Jf , evaluated at the solution t*, is roughly 0.89 for 
the first case (£2 = 0.16) and 1.02 for the second case 
(t 2 = 0.24). This explains why the first iteration con¬ 
verges, but not the second. Some further investigation 
reveals that this is due in large part to the discontinuities 
d fi 

in —. In our experiment the reference materials were 
at 

air, fat, soft tissue and bone, with tabulated LAC val¬ 
ues at the reference energy of 70 keV equal to 0, 0.1782, 
0.2033 and 0.4948 cm -1 , respectively. Thus the partial 
d (i 

derivative — has a larger value for t 2 = 0.24 (which lies 


Polyenergetic case Polyenergetic case 



FIG. 5. pSART iterations for t 2 = 0.203 (left) and t 2 = 0.204 
(right), narrowly centred on the solution. The first set of 
iterations converges to solution while the second converges to 
a two-cycle. 


between soft tissue and bone) than for t 2 = 0.16 (which 
lies between fat and air). In Fig. 5 we show the result of 
two more experiments where t 2 was set to values of 0.203 
and 0.204, which lie on either side of the tabulated value 
for soft tissue, where the derivative is discontinuous. The 
pSART iteration converges to the true solution in the first 
case but not in the second case, where it reaches a two- 
cycle between two points lying close to the true solution. 
Direct computation of the spectral radius confirms that 
it is equal to 0.87 in the first case and 1.16 in the second. 

Fig. 6 gives a convergence map for this test case as a 
function of t\ and t 2 . The effect of the discontinuities 
d fi 

in — is clearly visible in the discontinuities in p{Jf) 

that occur at the values of 0.1782 - the LAC value of 
fat at 70 keV - and 0.2033, the LAC value of soft tissue 
at 70 keV. It is apparent that the iteration transitions 
between convergent and non-convergent states at values 
of (ti,t 2 ) that do not he along these discontinuities as 
well. We note that the figure is not symmetric along 
the line t\ = t 2 - for example, all of the test cases that 
have been considered in this section would converge if 
the values of <1 and t 2 were interchanged. This figure 
is specific to the system matrix that arises from the ray 
paths illustrated in Fig. 1, and would be different for 
other paths. 


III. NUMERICAL EXPERIMENTS 

We have shown in the previous section that the pSART 
iteration is not guaranteed to converge in general, despite 
the success of the algorithm demonstrated in Ref. 18. 
One possible explanation is that the system matrices that 
arise in CT imaging are typically quite sparse and struc¬ 
tured, compared to the 2x2 matrix that was considered 
in the previous section. Thus, it is worth investigating 
whether the spectral radius of the Jacobian matrix of 
pSART is likely to exceed one for more realistic CT sys¬ 
tem matrices. 

In the following numerical experiment we consider the 
problem of reconstructing an TV x N pixel image for N = 
100, 200, 400 and 800. We simulate parallel beam data 
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P0 F ) Convergence Map 



FIG. 6. Spectral radius and convergence maps for the 2x2 
polyenergetic experiment. Left figure shows the spectral ra¬ 
dius p(Jf) as a function of t\ and t2‘, right figure is the equiv¬ 
alent binary map indicating whether the iteration converges 
(blue) or not (red). White crosses indicate the cases shown 
in Figs. 4 and 5. 



FIG. 7. FORBILD phantom for the case N = 800. Color 
window is restricted to [0.195, 0.215] to show low-contrast fea¬ 
tures; bony structures (white) have a LAC of 0.4948 at the 
reference energy of 70 keV. 

acquired at m equally spaced views over 180°, with m = 
180, 360, 720 and 1440, respectively. Forward projection 
(multiplication by A) is implemented using the radon 
command in Matlab, while backward projection (multi¬ 
plication by A T ) uses iradon with no filtering. Since 
Jf{ t*) depends on the object t* that we wish to recon¬ 
struct, we must consider a specific object to analyze the 
convergence of pSART. We use an N x N slice of the 
FORBILD numerical head phantom [22], which consists 
of bone and soft tissue, and includes some low-contrast 
features. An image of the phantom for the case N = 800 
is shown in Fig. 7. 

With these elements in place, we can approximate the 
spectral radius of the matrices associated with the SART 
iteration (11) and the Jacobian matrix of the pSART 
iteration (14, 15). To approximate the spectral radius, we 
use the power iteration (see for instance Ref. 23), which 
provides iterative estimates of the largest eigenvalue of a 
matrix M and the associated eigenvector, denoted by A ^ 
and x.( k \ respectively. The iterations were started from 
a random unit vector and run until the the largest 
element of the residual, ||Mx^ — A^x^Hoc, was less 
than 10” 4 . 

Table I shows the computed approximations to p(T) 


TABLE I. Results of the power iteration for the matrix T 
used for the SART iteration, and the Jacobian matrix of the 
pSART iteration, for different phantom sizes. N is the di¬ 
mension the image to be reconstructed and m is the total 
number of views. p(T) is the spectral radius of the matrix 
T, and p(Jf( t*)) is the spectral radius of the Jacobian of 
pSART iteration, evaluated for the N x N FORBILD phan¬ 
tom. The number of iterations run for each power iteration 
are displayed below the estimate of the spectral radius. 


N 

m 

p(T) p(J F ( t*)) 

100 

180 

0.999957 

0.999960 



2535 its 

2802 its 

200 

360 

0.999957 

0.999958 



2899 its 

2880 its 

400 

720 

0.999952 

0.999954 



3475 its 

3544 its 

800 

1440 

0.999933 

0.999935 



3544 its 

3595 its 


and p(Jf( t*)) for the different test cases, along with the 
number of power iterations required to obtain the esti¬ 
mate. It is apparent that there is virtually no difference 
between the spectral radius of the SART iteration matrix, 
T, and that of the the pSART Jacobian matrix, Jp(t*), 
in any of the studied cases. In no instances did the com¬ 
puted estimate of ever exceed one, which would 

cause the iteration to diverge in the neighbourhood of the 
solution. We conclude that for more realistic CT imag¬ 
ing scenarios, the pSART iteration is, at the very least, 
likely to exhibit local convergence in the neighbourhood 
of the solution, with a comparable rate of convergence 
to SART. We note, however, that the spectral radius of 
both the SART and pSART iterations is very close to one, 
indicating that the convergence will be slow in the neigh¬ 
bourhood of the solution. The convergence can likely be 
accelerated by using subsets of the projection data, as 
was proposed in the original SART algorithm [4], 

This analysis establishes only local convergence in the 
neighbourhood of the solution. Global convergence (i.e. 
from an arbitrary initial estimate) is more difficult to 
establish in general, and it is not obvious whether there 
exist conditions on the system matrix A , beam spectrum, 
choice of reference energy, etc. to guarantee global con¬ 
vergence of pSART. This is fairly typical of nonlinear 
fixed-point iterations, Newton’s method being a well- 
known example. The method does appear to be fairly 
robust with respect to the choice of starting point, how¬ 
ever, as the images in Ref. 18 were produced from an 
initial estimate consisting only of zeros; our own expe¬ 
rience with the algorithm confirms that this choice of 
initial estimate works well. 

Additionally, our analysis assumes the existence of an 
exact solution, whereas in practice the system is likely 
to be inconsistent due to factors such as noisy measure¬ 
ments and model mismatch. As with any reconstruction 














6 


algorithm, model mismatch will produce artifacts in the 
reconstructed image; Ref. 18 includes some experiments 
quantifying the effect of spectrum error. In the presence 
of noisy data, our experience indicates that pSART ex¬ 
hibits the “semi-convergence” behaviour typical of other 
iterative methods (see e.g. Ref. 24, p.89); namely, that 
the algorithm initially converges towards the solution, 
but the image eventually deterioriates with further iter¬ 
ations due to the effects of noise. As noted in Ref. 18, 
this problem could potentially be addressed with the use 
of statistical modeling or edge-preserving regularization. 


CONCLUSIONS 

In this paper we have analyzed the convergence of a 
recently proposed polyenergetic SART (pSART) algo¬ 
rithm. We show that the spectral radius of the Jacobian 
of the nonlinear pSART iteration may be larger than one 
in some cases. Thus the method is not mathematically 
guaranteed to converge to a solution of the nonlinear sys¬ 
tem of polyenergetic equations, in general. For system 
matrices of the type encountered in CT imaging, how¬ 
ever, our empirical results indicate that the spectral ra¬ 
dius of the Jacobian matrix, evaluated for a prototypical 
head phantom, is essentially the same as the spectral ra¬ 
dius of the matrix corresponding to the convergent SART 
iteration. Thus in practice it seems that the method is 
likely to converge at roughly the same rate as SART. 
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Lemma A.l: Let W = DA T MA. Then, p(W) < 1, 
with equality if all elements of A are positive. 

Proof: A direct calculation gives the (i,j )th element of 
W as 


Wij 


1 

Ji 


E 


k =1 


Q j kiQ j kj 
7 k 


It follows that the sum of the absolute values in row i of 
W is: 


rti 

1 x > &kiQ'kj 


Ik 

3 — 1 j—l k—1 


1 m | | n 

dEyEw 

Pl k=1 lk i=l 


_ 1 \ a ki\ _ 

7fc 

1 m 


-7k 


k=1 


-S* 

= 1 


Since the max norm of a matrix, || • Hoc, is equal to the 
maximum row sum, and the spectral radius of a matrix 
cannot be greater than the max norm, it must be true 
that 


p { w ) < Halloo < i. 
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APPENDIX 

We provide a short proof of convergence for SART sub¬ 
ject to a condition on the rank of A. Convergence of 
SART has been established previously in Refs. 6 and 
7. This proof is somewhat shorter and is intended to 
complement the analysis we have presented for pSART. 

Let T = I — DA T MA, where A is the m x n system 
matrix and D and M are defined in (4). We assume that 
A has rank n, meaning that there are as many linearly 
independent equations as there are unknowns. This im¬ 
plies that m > n, i.e. that there are at least as many 
measurements as there are unknowns. The SART itera¬ 
tion then has the form 

x (fc+i) = Tx (U + c 

We first prove two lemmas. 


When all elements of A are positive (which is the case in 
tomographic applications), every row of W sums exactly 
to 1, and so A = 1 is an eigenvalue of W (with the asso¬ 
ciated eigenvector consisting of all ones), and p{W) = 1. 

On its own, this lemma only tells us that the eigenval¬ 
ues of W have magnitude less than 1. We also need the 
following result: 

Lemma A.2: All eigenvalues of W are positive real 
numbers. 

Proof: This Lemma is an application of Theorem 7.6.3 
from Ref. 25. We first remind the reader of two defini¬ 
tions: 

1. Two nxn matrices A and B are similar if there ex¬ 
ists an invertible matrix P such that A = PBP -1 . 
If A and B are similar, then they have exactly the 
same eigenvalues. 

2. Two real nxn matrices A and B are congruent if 
there exists an invertible matrix P such that A = 
PBP T . If A and B are symmetric, then Sylvester’s 
law of inertia states that they have the same inertia, 
meaning the same number of positive, negative, and 
zero eigenvalues. (Recall that the eigenvalues of a 
symmetric matrix are always real numbers). 
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Now, let V = A T MA. Then, V is a symmetric pos¬ 
itive semidefinite n x n matrix, since it can be writ¬ 
ten as the product of a matrix and its transpose: V = 
(M 5 A) t (M 5 A). Furthermore, since A has rank n and 
M 5 is a diagonal matrix, M 2 A has rank n and so does 
V. It follows that V is invertible, and so zero is not an 
eigenvalue of V. Thus all eigenvalues of V are positive. 
We then have the following two results: 

1. W = DV is similar to D^VD?, since D^VD^ = 

D-^DVD^. 

2. D^VD^ is congruent to V. 


The first result implies that the eigenvalues of W are the 
same as the eigenvalues of D 5 VD 5 , while the second im¬ 
plies that the eigenvalues of this matrix must have the 
same signs as the eigenvalues of V. So, since all eigen¬ 
values of V are positive, all eigenvalues of W must be 
positive as well. 

Theorem A.3: The matrix T satisfies p(T) < 1, and 
hence the SART iteration converges. 

Proof: Lemmas A.l and A.2 prove that any eigenvalues 
A of IT satisfy 0 < A < 1. Thus the eigenvalues of T = 
I — W satisfy 0 < A < 1, and so p{T) < 1. 
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